FUT2_Baseline
Structure and composition
Yiming Wang
30/05/2024
Reference sites
Packages
pacman::p_load(
dplyr,
tidyverse,
readxl,
vegan)
# setwd("D:/PhD/01_Data/01_Baseline/Model_1") # Windows
# setwd("/Volumes/Yiming_Wang/PhD/01_Data/01_Baseline/Model_1") # Mac
# setwd("/Users/godric/Desktop/Fut2 review comments/Revision") #Yiming's macbookpro
setwd("/Volumes/Yiming_Wang/PhD/15_Manuscript/Submission/07_ISME/Review comments/Revision data/26 April 2024")#Yiming's imac pro
df_all <- read_excel("Baseline_S33_taxa_formatted.xlsx")
# Subgroup by variables
#30436d WT Male
#9ea5c2 WT Female
#654e3a KO Male
#bcaba6 KO Female
#3C5488B2 WT
#7E6148B2 KO
#D98594FF Female
#94C5CCFF MaleAlpha diversity_visualization and analysis
Plot-all mice
# load package
pacman::p_load(
ggplot2,ggsci, ggthemes,
)
set.seed(123)
# Level orders
level_order_genotype <- c("WT","KO")
level_order_sex <- c("Male","Female")
level_order_genotypesex <- c("WT Male","WT Female","KO Male","KO Female")
# My color panel
mypal = pal_npg("nrc", alpha = 0.7)(9)
mypal## [1] "#E64B35B2" "#4DBBD5B2" "#00A087B2" "#3C5488B2" "#F39B7FB2" "#8491B4B2"
## [7] "#91D1C2B2" "#DC0000B2" "#7E6148B2"library("scales")
show_col(mypal) # #3C5488B2(WT) and #E64B35B2 (KO)
# Genotype
## FDP Diversity
Genotype_FDP <- ggplot (data=df_all, mapping = aes(y=Diversity_FPD, x=factor(Genotype, level=level_order_genotype))) +
geom_boxplot(outlier.shape = NA)+
theme_bw() +
geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
position=position_jitterdodge(jitter.width = 0.5,
jitter.height = 0,
dodge.width = 0.8))+
labs(x="", y="Faith's phylogenetic diversity", caption ="", # x-,y-axis name
title = "")+
theme(legend.position = "none", # right,left, bottom
legend.title = element_blank(), # can remove this
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 5)),
axis.title.y=element_text(margin = margin(r = 5), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_y_continuous(limits = c(4,24),breaks = seq(4,24,5))+
scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) +
annotate("text",x=1.5,y=24,
label="ns",hjust=0.5, size=5.5)
Genotype_FDP_V <- Genotype_FDP + geom_violin(alpha=0, width=0.9,
position=position_dodge(width = 0.8),
size=0.75)
## Shannon Diversity
Genotype_shannon <- ggplot (data=df_all, mapping = aes(y=Diversity_Shannon, x=factor(Genotype, level=level_order_genotype))) +
geom_boxplot(outlier.shape = NA)+
theme_bw() +
geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
position=position_jitterdodge(jitter.width = 0.5,
jitter.height = 0,
dodge.width = 0.8))+
labs(x="", y="Shannon diversity", caption ="", # x-,y-axis name
title = "")+
theme(legend.position = "none", # right,left, bottom
legend.title = element_blank(), # can remove this
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 5)),
axis.title.y=element_text(margin = margin(r = 5), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_y_continuous(limits = c(5,8),breaks = seq(5,8,1))+
scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) +
annotate("text",x=1.5,y=8,
label="ns",hjust=0.5, size=5.5)
Genotype_shannon_V <- Genotype_shannon + geom_violin(alpha=0, width=0.9,
position=position_dodge(width = 0.8),
size=0.75)
## Richness (Chao)
Genotype_chao <- ggplot (data=df_all, mapping = aes(y=Richness_Chaos, x=factor(Genotype, level=level_order_genotype))) +
geom_boxplot(outlier.shape = NA)+
theme_bw() +
geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
position=position_jitterdodge(jitter.width = 0.5,
jitter.height = 0,
dodge.width = 0.8))+
labs(x="", y="Chao richness", caption ="", # x-,y-axis name
title = "")+
theme(legend.position = "none", # right,left, bottom
legend.title = element_blank(), # can remove this
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 5)),
axis.title.y=element_text(margin = margin(r = 5), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_y_continuous(limits = c(50,250),breaks = seq(50,250,50))+
scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) +
annotate("text",x=1.5,y=250,
label="ns",hjust=0.5, size=5.5)
Genotype_chao_V <- Genotype_chao + geom_violin(alpha=0, width=0.9,
position=position_dodge(width = 0.8),
size=0.75)
## Pielou evenness
Genotype_pielou <- ggplot (data=df_all, mapping = aes(y=Evenness_Pielou, x=factor(Genotype, level=level_order_genotype))) +
geom_boxplot(outlier.shape = NA)+
theme_bw() +
geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
position=position_jitterdodge(jitter.width = 0.5,
jitter.height = 0,
dodge.width = 0.8))+
labs(x="", y="Pielou evenness", caption ="", # x-,y-axis name
title = "")+
theme(legend.position = "none", # right,left, bottom
legend.title = element_blank(), # can remove this
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 5)),
axis.title.y=element_text(margin = margin(r = 5), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_y_continuous(limits = c(0.8,1.0),breaks = seq(0.8,1.0,0.05))+
scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) +
annotate("text",x=1.5,y=1.0,
label="ns",hjust=0.5, size=5.5)
Genotype_pielou_V <- Genotype_pielou + geom_violin(alpha=0, width=0.9,
position=position_dodge(width = 0.8),
size=0.75)
# Print
Genotype_FDP
Genotype_shannon
Genotype_chao
Genotype_pielou
Plot-male mice only
# Level orders
level_order_genotype <- c("WT","KO")
level_order_sex <- c("Male","Female")
level_order_genotypesex <- c("WT Male","WT Female","KO Male","KO Female")
# My color panel
mypal = pal_npg("nrc", alpha = 0.7)(9)
mypal## [1] "#E64B35B2" "#4DBBD5B2" "#00A087B2" "#3C5488B2" "#F39B7FB2" "#8491B4B2"
## [7] "#91D1C2B2" "#DC0000B2" "#7E6148B2"library("scales")
show_col(mypal) # #3C5488B2(WT) and #E64B35B2 (KO)
df_male <- df_all %>%
filter(Sex == "Male")
## FDP Diversity
Genotype_FDP_male <- ggplot (data=df_male, mapping = aes(y=Diversity_FPD, x=factor(Genotype, level=level_order_genotype))) +
geom_boxplot(outlier.shape = NA)+
theme_bw() +
geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
position=position_jitterdodge(jitter.width = 0.5,
jitter.height = 0,
dodge.width = 0.8))+
labs(x="", y="Faith's phylogenetic diversity", caption ="", # x-,y-axis name
title = "")+
theme(legend.position = "none", # right,left, bottom
legend.title = element_blank(), # can remove this
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 5)),
axis.title.y=element_text(margin = margin(r = 5), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_y_continuous(limits = c(4,24),breaks = seq(4,24,5))+
scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) +
annotate("text",x=1.5,y=24,
label="ns",hjust=0.5, size=5.5)
Genotype_FDP_V <- Genotype_FDP + geom_violin(alpha=0, width=0.9,
position=position_dodge(width = 0.8),
size=0.75)
## Shannon Diversity
Genotype_shannon_male <- ggplot (data=df_male, mapping = aes(y=Diversity_Shannon, x=factor(Genotype, level=level_order_genotype))) +
geom_boxplot(outlier.shape = NA)+
theme_bw() +
geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
position=position_jitterdodge(jitter.width = 0.5,
jitter.height = 0,
dodge.width = 0.8))+
labs(x="", y="Shannon diversity", caption ="", # x-,y-axis name
title = "")+
theme(legend.position = "none", # right,left, bottom
legend.title = element_blank(), # can remove this
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 5)),
axis.title.y=element_text(margin = margin(r = 5), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_y_continuous(limits = c(5,8),breaks = seq(5,8,1))+
scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) +
annotate("text",x=1.5,y=8,
label="ns",hjust=0.5, size=5.5)
Genotype_shannon_V <- Genotype_shannon + geom_violin(alpha=0, width=0.9,
position=position_dodge(width = 0.8),
size=0.75)
## Richness (Chao)
Genotype_chao_male <- ggplot (data=df_male, mapping = aes(y=Richness_Chaos, x=factor(Genotype, level=level_order_genotype))) +
geom_boxplot(outlier.shape = NA)+
theme_bw() +
geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
position=position_jitterdodge(jitter.width = 0.5,
jitter.height = 0,
dodge.width = 0.8))+
labs(x="", y="Chao richness", caption ="", # x-,y-axis name
title = "")+
theme(legend.position = "none", # right,left, bottom
legend.title = element_blank(), # can remove this
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 5)),
axis.title.y=element_text(margin = margin(r = 5), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_y_continuous(limits = c(50,250),breaks = seq(50,250,50))+
scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) +
annotate("text",x=1.5,y=250,
label="ns",hjust=0.5, size=5.5)
Genotype_chao_V <- Genotype_chao + geom_violin(alpha=0, width=0.9,
position=position_dodge(width = 0.8),
size=0.75)
## Pielou evenness
Genotype_pielou_male <- ggplot (data=df_male, mapping = aes(y=Evenness_Pielou, x=factor(Genotype, level=level_order_genotype))) +
geom_boxplot(outlier.shape = NA)+
theme_bw() +
geom_jitter(aes(colour=Genotype),size=2,alpha=0.9,
position=position_jitterdodge(jitter.width = 0.5,
jitter.height = 0,
dodge.width = 0.8))+
labs(x="", y="Pielou evenness", caption ="", # x-,y-axis name
title = "")+
theme(legend.position = "none", # right,left, bottom
legend.title = element_blank(), # can remove this
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 5)),
axis.title.y=element_text(margin = margin(r = 5), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_y_continuous(limits = c(0.8,1.0),breaks = seq(0.8,1.0,0.05))+
scale_colour_manual(values=c("#7E6148B2","#3C5488B2")) +
annotate("text",x=1.5,y=1.0,
label="ns",hjust=0.5, size=5.5)
Genotype_pielou_V <- Genotype_pielou + geom_violin(alpha=0, width=0.9,
position=position_dodge(width = 0.8),
size=0.75)
# Print
Genotype_FDP_male
Genotype_shannon_male
Genotype_chao_male
Genotype_pielou_male
Beta diversity_results and visualisation
Datashape
# places that are not satisfied
## 1) Permanova results table
## 2) Figure A and B did not align because of legend (Change legend size, position to bottom, change P value position)
#library packages
pacman::p_load(
data.table,
vegan,
knitr,
EcolUtils
)
# Shape data
# Level orders
df_all$Genotype <- factor (df_all$Genotype, levels = c("WT", "KO"))
df_all$Sex <- factor(df_all$Sex, levels = c("Male", "Female"))
df_all$GenotypeSex <- factor(df_all$GenotypeSex, levels = c("WT Male", "WT Female", "KO Male","KO Female"))
# subset dateframe
df_abundance <- subset(df_all, select = -c(1:17))
df_metadata <- subset(df_all, select = c(1:17))Bray-curtis
set.seed(123)
# Check the stress and choose the optimal number of dimensions
dist_bray <- vegdist(df_abundance, method = "bray")
# NMDS.scree(dist)
## For a good representation of your data, the acceptable stress value: <0.2
## Considering our stress value is 0.13, which is good enough for representation of microbiome composition, so k=2 is chosen for NMDS plot
## The stress value does not affect PERMANOVA results, only tells you whether 2D (k=2) is good for composition demonstration
# NMDS ordination_bray curtis
NMDS <- metaMDS(df_abundance, distance = "bray", k = 2)## Run 0 stress 0.1413421
## Run 1 stress 0.1415037
## ... Procrustes: rmse 0.004186663 max resid 0.02572238
## Run 2 stress 0.1416862
## ... Procrustes: rmse 0.007319768 max resid 0.04434329
## Run 3 stress 0.1413421
## ... Procrustes: rmse 2.594225e-05 max resid 0.0001363623
## ... Similar to previous best
## Run 4 stress 0.1413421
## ... Procrustes: rmse 5.953767e-06 max resid 3.293415e-05
## ... Similar to previous best
## Run 5 stress 0.1413421
## ... Procrustes: rmse 3.72446e-05 max resid 0.0002456958
## ... Similar to previous best
## Run 6 stress 0.14154
## ... Procrustes: rmse 0.005891884 max resid 0.04412429
## Run 7 stress 0.1415037
## ... Procrustes: rmse 0.004208291 max resid 0.02583719
## Run 8 stress 0.1416861
## ... Procrustes: rmse 0.007331071 max resid 0.04434624
## Run 9 stress 0.1844246
## Run 10 stress 0.1615297
## Run 11 stress 0.161271
## Run 12 stress 0.1413422
## ... Procrustes: rmse 8.19954e-05 max resid 0.0004571962
## ... Similar to previous best
## Run 13 stress 0.1413421
## ... New best solution
## ... Procrustes: rmse 4.073193e-06 max resid 2.546618e-05
## ... Similar to previous best
## Run 14 stress 0.1612711
## Run 15 stress 0.1413421
## ... Procrustes: rmse 1.713322e-05 max resid 0.0001087092
## ... Similar to previous best
## Run 16 stress 0.1413421
## ... Procrustes: rmse 2.929669e-05 max resid 0.0001228787
## ... Similar to previous best
## Run 17 stress 0.1416861
## ... Procrustes: rmse 0.007329494 max resid 0.04434847
## Run 18 stress 0.1413421
## ... Procrustes: rmse 1.40533e-05 max resid 9.236742e-05
## ... Similar to previous best
## Run 19 stress 0.1413423
## ... Procrustes: rmse 3.979513e-05 max resid 0.000235021
## ... Similar to previous best
## Run 20 stress 0.1415037
## ... Procrustes: rmse 0.004176246 max resid 0.02566424
## *** Solution reachedNMDS$stress## [1] 0.1413421# PERMANOVA and Pair-wise PERMANOVA
set.seed(123)
## Genotype and sex are crossed variable, cage are nested in genotype and sex: https://stats.stackexchange.com/questions/350462/can-you-perform-a-permanova-analysis-on-nested-data
adonis2(df_abundance ~ (Genotype+Sex)/Cage, data=df_metadata, permutations=9999, method="bray")## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = df_abundance ~ (Genotype + Sex)/Cage, data = df_metadata, permutations = 9999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Genotype 1 0.1057 0.02780 2.7448 0.0280 *
## Sex 1 0.0824 0.02166 2.1391 0.0721 .
## Genotype:Sex:Cage 16 1.8816 0.49481 3.0537 0.0001 ***
## Residual 45 1.7330 0.45573
## Total 63 3.8027 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Genotype x Sex effect
adonis2(df_abundance ~ GenotypeSex, data=df_metadata, permutations=9999, method="bray") ## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 9999
##
## adonis2(formula = df_abundance ~ GenotypeSex, data = df_metadata, permutations = 9999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## GenotypeSex 3 0.4048 0.10645 2.3827 0.0068 **
## Residual 60 3.3979 0.89355
## Total 63 3.8027 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1df_permanova_pairwise <- adonis.pair(vegdist(df_all[,18:128], method = "bray"), df_all$GenotypeSex, nper = 9999, corr.method = "fdr")
kable(df_permanova_pairwise)| combination | SumsOfSqs | MeanSqs | F.Model | R2 | P.value | P.value.corrected |
|---|---|---|---|---|---|---|
| WT Male <-> WT Female | 0.0700924 | 0.0700924 | 1.2445858 | 0.0398336 | 0.2827 | 0.38328 |
| WT Male <-> KO Male | 0.2635370 | 0.2635370 | 4.2755571 | 0.1247407 | 0.0041 | 0.02130 |
| WT Male <-> KO Female | 0.0521060 | 0.0521060 | 0.9933835 | 0.0320515 | 0.4059 | 0.40590 |
| WT Female <-> KO Male | 0.1359739 | 0.1359739 | 2.2360774 | 0.0693657 | 0.0607 | 0.12140 |
| WT Female <-> KO Female | 0.0588837 | 0.0588837 | 1.1406231 | 0.0366281 | 0.3194 | 0.38328 |
| KO Male <-> KO Female | 0.2290033 | 0.2290033 | 4.0215307 | 0.1182055 | 0.0071 | 0.02130 |
# Assessment of multivariate homogeneity of groups dispersions
dispersion_bray <- betadisper(dist_bray, group = df_metadata$GenotypeSex)
permutest(dispersion_bray)##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 3 0.00865 0.0028848 0.4361 999 0.725
## Residuals 60 0.39688 0.0066146# Figures
#library packages
pacman::p_load(
ggthemes,
ggsci
)
## Set seed
set.seed(123)
# options( scipen = 999 )
# Genotype
datascores <- as.data.frame(scores(NMDS))# $sites
scores <- cbind(as.data.frame(datascores), Genotype = df_metadata$Genotype)
centroids <- aggregate(cbind(NMDS1, NMDS2) ~ Genotype, data = scores, FUN = mean)
seg <- merge(scores, setNames(centroids, c('Genotype','oNMDS1','oNMDS2')),
by = c('Genotype'), sort = FALSE)
composition_Genotype <- ggplot(scores, aes(x = NMDS1, y = NMDS2, color = Genotype)) +
geom_segment(data = seg,
mapping = aes(xend = oNMDS1, yend = oNMDS2),
size=0.25) +
geom_point(data = centroids, size = 4) +
geom_point()+
coord_fixed()+
theme_bw() +
labs(title="Genotype (Distance matrix: Bray-Curtis)")+
theme(legend.position = "right",
legend.title= element_text(size=10,face="plain"),
legend.text = element_text(size=10),
legend.key = element_rect(fill = "transparent", colour = "black"),
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 5)),
axis.title.y=element_text(margin = margin(r = 5), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_x_continuous(limits = c(-0.5,0.75),breaks = seq(-0.5,0.75,0.25))+
scale_y_continuous(limits = c(-0.5,0.75),breaks = seq(-0.5,0.75,0.25))+
scale_color_manual(values=c("#3C5488B2","#7E6148B2"))+
scale_fill_manual(values=c("#3C5488B2","#7E6148B2"))+
annotate("text",x=-0.5,y=0.75, #0.4 if legend position is right
label="Stress = 0.14",hjust=0, size=3)+ #2.2 if legend position is right
annotate("text",x=-0.5,y=0.69,
label="Pseudo-F = 2.74",hjust=0, size=3)+
annotate("text",x=-0.5,y=0.63,
label="p = 0.028",hjust=0, size=3)
# Genotype x Sex
datascores <- as.data.frame(scores(NMDS)) #$sites
scores <- cbind(as.data.frame(datascores), GenotypeSex = df_metadata$GenotypeSex) %>%
separate(GenotypeSex,c("Genotype","Sex"),sep="_",remove=FALSE)
centroids <- aggregate(cbind(NMDS1, NMDS2) ~ GenotypeSex, data = scores, FUN = mean) %>%
separate(GenotypeSex,c("Genotype","Sex"),sep="_",remove=FALSE)
seg <- merge(scores, setNames(centroids, c('GenotypeSex','Genotype','Sex','oNMDS1','oNMDS2')),
by = c('GenotypeSex','Genotype','Sex'), sort = FALSE)
composition_GenotypeSex_final <- ggplot(scores, aes(x = NMDS1, y = NMDS2,color= GenotypeSex,fill=GenotypeSex)) + #, colour =GenotypeSex
geom_segment(data = seg,
mapping = aes(xend = oNMDS1, yend = oNMDS2),
size = 0.25) +
geom_point(data = centroids, size = 4) +
geom_point()+ #geo,_point()
coord_fixed()+
theme_bw() +
labs(title="")+
theme(legend.position = "bottom",
legend.title= element_text(size=10, face="plain"),
legend.text = element_text(size=14),
legend.key = element_rect(fill = "transparent", colour = "black"),
axis.title = element_text(face="plain", size = 10),
axis.title.x=element_text(margin = margin(t = 9), size=16),
axis.title.y=element_text(margin = margin(r = 3), size=16),
axis.text.x=element_text(colour="black", face="plain", size=14),
axis.text.y=element_text(colour="black", face="plain", size=14),
plot.title = element_text(size=10, hjust=0.5, face="plain"),
panel.border = element_rect(color="black", #element_rect is often used for backgrounds and borders
fill = NA,
size = 1),
legend.background = element_rect(color="transparent"),
panel.grid = element_blank(),
aspect.ratio = 1,
plot.background = element_rect(fill="transparent", colour = "transparent"))+
scale_x_continuous(limits = c(-0.5,0.75),breaks = seq(-0.5,0.75,0.25))+
scale_y_continuous(limits = c(-0.5,0.75),breaks = seq(-0.5,0.75,0.25))+
scale_color_manual(values=c("#30436d","#9ea5c2","#654e3a","#bcaba6"))+
scale_fill_manual(values=c("#30436d","#9ea5c2","#654e3a","#bcaba6"))+
scale_shape_manual(values=c(24,21,24,21))+
annotate("text",x=-0.5,y=-0.5,
label="Stress = 0.14",hjust=0, size=4)+
annotate("text",x=-0.5,y=0.75,
label=expression(paste("Male (WT vs KO): ",italic("P") == 0.021)),hjust=0, size=4)+
annotate("text",x=-0.5,y=0.65,
label=expression(paste("Female (WT vs KO): ",italic("P") > 0.05)),hjust=0, size=4)
# Print
composition_GenotypeSex_final
ANOISM: Variance between groups and within groups
- ANOSIM analysis to test whether the variance between groups is larger than within groups (if it is larger, then P values is less than 0.05)
- Anoism analysis (R value from -1 to 1, R > 0 indicates between-group difference > within group difference; R < 0 indicates within group > between group difference; P value indicates the statistical significance)
# Bray-curtis
dist_bray <- vegdist(df_abundance, method = "bray")
## Genotype
anosim.genotype_bray <- with(df_metadata, anosim(dist_bray, Genotype))
summary(anosim.genotype_bray)##
## Call:
## anosim(x = dist_bray, grouping = Genotype)
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.02295
## Significance: 0.111
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0260 0.0407 0.0556 0.0681
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 2 515.75 1035.0 1513.50 2016 1024
## WT 3 506.75 944.0 1453.25 2015 496
## KO 1 481.00 1028.5 1596.25 2011 496plot(anosim.genotype_bray)
## Genotype x Sex
anosim.GenotypeSex_bray <- with(df_metadata, anosim(dist_bray, GenotypeSex))
summary(anosim.GenotypeSex_bray)##
## Call:
## anosim(x = dist_bray, grouping = GenotypeSex)
## Dissimilarity: bray
##
## ANOSIM statistic R: 0.0745
## Significance: 0.014
##
## Permutation: free
## Number of permutations: 999
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.0333 0.0449 0.0570 0.0770
##
## Dissimilarity ranks between and within classes:
## 0% 25% 50% 75% 100% N
## Between 2 532.75 1031.5 1521.25 2016 1536
## WT Male 3 349.00 784.5 1512.25 2010 120
## WT Female 7 619.75 969.5 1431.00 1934 120
## KO Male 5 457.50 1154.0 1704.25 1979 120
## KO Female 1 395.00 816.5 1192.25 1965 120plot(anosim.GenotypeSex_bray)